###########################################################
### Disorganized Political Violence:                    ###
### A Demonstration Case of Temperature and Insurgency  ###
### Replication Code:                                   ###
### Temperature Response Function Results:              ### 
### Country-Wide Robustness Checks                      ###
### Andrew Shaver, Alex Bollfrass                       ###
### Date: 07/04/2022                                    ###
###########################################################

# Note: these results were generated using:
# 1) MacBook Pro (16-inch, 2021); Chip: Apple M1 Max; Memory: 64 GB; running
# macOS Monterey.
# 2) R 4.1.2 GUI 1.77 High Sierra build (8007)
# 3) RStudio 2022.07.1+554 "Spotted Wakerobin" Release (7872775ebddc40635780ca1ed238934c3345c5de, 2022-07-22) for macOS
# Mozilla/5.0 (Macintosh; Intel Mac OS X 12_1_0) AppleWebKit/537.36 (KHTML, like Gecko) QtWebEngine/5.12.10 Chrome/69.0.3497.128 Safari/537.36
# (Please see below for specific package versions.)

# Load required packages:
library(broom)
library(dplyr)
library(stargazer)
library(sandwich)
library(lfe)
library(ggplot2)

packageVersion("broom")
# ‘0.7.12’
packageVersion("dplyr")
# ‘1.0.9’
packageVersion("stargazer")
# ‘5.2.3’
packageVersion("sandwich")
# ‘3.0.1’
packageVersion("lfe")
# ‘2.8.7.1’
packageVersion("ggplot2")
# ‘3.3.6’

# Load datasets:
panel_af <- readRDS("~/Library/CloudStorage/Box-Box/TurningUpTheHeat1/IO_Replication_Materials/Processed Data/af_panel.Rds")
panel_iq <- readRDS("~/Library/CloudStorage/Box-Box/TurningUpTheHeat1/IO_Replication_Materials/Processed Data/iq_panel.Rds")

all_observed_temps <- c(panel_af$temperature, panel_iq$temperature)
all_observed_temps <- as.data.frame(all_observed_temps)
colnames(all_observed_temps) <- "temperature"

panel_af$temperature_cubed <- panel_af$temperature^{3}
panel_iq$temperature_cubed <- panel_iq$temperature^{3}

# Plot predicted values across various specifications:
# Find minimum and maximum temperatures:
temperature_range <- round(min(panel_af$temperature, panel_iq$temperature):round(max(panel_af$temperature, panel_iq$temperature)))
quantile(temperature_range, c(0.025, 0.975))
temperature_range <- temperature_range[temperature_range %in% 
                                         c(round(as.numeric(quantile(temperature_range, c(0.025, 0.975))[1])):round(as.numeric(quantile(temperature_range, c(0.025, 0.975))[2])))]

AF.w.g.t.l <- felm(df ~ temperature  + precipitation + idf + ied +
                     temperature_1 + temperature_2 + temperature_3 + temperature_4 + temperature_5 + temperature_6 + temperature_7 +
                     precipitation_1 + precipitation_2 + precipitation_3 + precipitation_4 + precipitation_5 + precipitation_6 + precipitation_7 +
                     df_1 + df_2 + df_3 + df_4 + df_5 + df_6 + df_7 +
                     idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                     ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                     week + DISTID | 0 | DISTID, data = panel_af)
IQ.w.g.t.l <- felm(unconstrained1 ~ temperature + precipitation + constrained1 + ied +
                     temperature_1 + temperature_2 + temperature_3 + temperature_4 + temperature_5 + temperature_6 + temperature_7 +
                     precipitation_1 + precipitation_2 + precipitation_3 + precipitation_4 + precipitation_5 + precipitation_6 + precipitation_7 +
                     unconstrained1_1 + unconstrained1_2 + unconstrained1_3 + unconstrained1_4 + unconstrained1_5 + unconstrained1_6 + unconstrained1_7 +
                     constrained1_1 + constrained1_2 + constrained1_3 + constrained1_4 + constrained1_5 + constrained1_6 + constrained1_7 +
                     ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                     week + ADM3NAME | 0 | ADM3NAME, data = panel_iq)


AF.w.g.t.q <- felm(df ~ temperature + temperature_squared  + precipitation + idf + ied +
                     temperature_1 + temperature_2 + temperature_3 + temperature_4 + temperature_5 + temperature_6 + temperature_7 +
                     precipitation_1 + precipitation_2 + precipitation_3 + precipitation_4 + precipitation_5 + precipitation_6 + precipitation_7 +
                     df_1 + df_2 + df_3 + df_4 + df_5 + df_6 + df_7 +
                     idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                     ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                     week + DISTID | 0 | DISTID, data = panel_af)
IQ.w.g.t.q <- felm(unconstrained1 ~ temperature + temperature_squared + precipitation + constrained1 + ied +
                     temperature_1 + temperature_2 + temperature_3 + temperature_4 + temperature_5 + temperature_6 + temperature_7 +
                     precipitation_1 + precipitation_2 + precipitation_3 + precipitation_4 + precipitation_5 + precipitation_6 + precipitation_7 +
                     unconstrained1_1 + unconstrained1_2 + unconstrained1_3 + unconstrained1_4 + unconstrained1_5 + unconstrained1_6 + unconstrained1_7 +
                     constrained1_1 + constrained1_2 + constrained1_3 + constrained1_4 + constrained1_5 + constrained1_6 + constrained1_7 +
                     ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                     week + ADM3NAME | 0 | ADM3NAME, data = panel_iq)

AF.w.g.t.c <- felm(df ~ temperature + temperature_squared + temperature_cubed  + precipitation + idf + ied +
                     temperature_1 + temperature_2 + temperature_3 + temperature_4 + temperature_5 + temperature_6 + temperature_7 +
                     precipitation_1 + precipitation_2 + precipitation_3 + precipitation_4 + precipitation_5 + precipitation_6 + precipitation_7 +
                     df_1 + df_2 + df_3 + df_4 + df_5 + df_6 + df_7 +
                     idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                     ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                     week + DISTID | 0 | DISTID, data = panel_af)
IQ.w.g.t.c <- felm(unconstrained1 ~ temperature + temperature_squared + temperature_cubed + precipitation + constrained1 + ied +
                     temperature_1 + temperature_2 + temperature_3 + temperature_4 + temperature_5 + temperature_6 + temperature_7 +
                     precipitation_1 + precipitation_2 + precipitation_3 + precipitation_4 + precipitation_5 + precipitation_6 + precipitation_7 +
                     unconstrained1_1 + unconstrained1_2 + unconstrained1_3 + unconstrained1_4 + unconstrained1_5 + unconstrained1_6 + unconstrained1_7 +
                     constrained1_1 + constrained1_2 + constrained1_3 + constrained1_4 + constrained1_5 + constrained1_6 + constrained1_7 +
                     ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                     week + ADM3NAME | 0 | ADM3NAME, data = panel_iq)

AF.m.g.t.l <- felm(df ~ temperature + precipitation + idf + ied +
                     temperature_1 + temperature_2 + temperature_3 + temperature_4 + temperature_5 + temperature_6 + temperature_7 +
                     precipitation_1 + precipitation_2 + precipitation_3 + precipitation_4 + precipitation_5 + precipitation_6 + precipitation_7 +
                     df_1 + df_2 + df_3 + df_4 + df_5 + df_6 + df_7 +
                     idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                     ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                     month + DISTID | 0 | DISTID, data = panel_af)
IQ.m.g.t.l <- felm(unconstrained1 ~ temperature + precipitation + constrained1 + ied +
                     temperature_1 + temperature_2 + temperature_3 + temperature_4 + temperature_5 + temperature_6 + temperature_7 +
                     precipitation_1 + precipitation_2 + precipitation_3 + precipitation_4 + precipitation_5 + precipitation_6 + precipitation_7 +
                     unconstrained1_1 + unconstrained1_2 + unconstrained1_3 + unconstrained1_4 + unconstrained1_5 + unconstrained1_6 + unconstrained1_7 +
                     constrained1_1 + constrained1_2 + constrained1_3 + constrained1_4 + constrained1_5 + constrained1_6 + constrained1_7 +
                     ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                     month + ADM3NAME | 0 | ADM3NAME, data = panel_iq)

AF.m.g.t.q <- felm(df ~ temperature + temperature_squared + precipitation + idf + ied +
                     temperature_1 + temperature_2 + temperature_3 + temperature_4 + temperature_5 + temperature_6 + temperature_7 +
                     precipitation_1 + precipitation_2 + precipitation_3 + precipitation_4 + precipitation_5 + precipitation_6 + precipitation_7 +
                     df_1 + df_2 + df_3 + df_4 + df_5 + df_6 + df_7 +
                     idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                     ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                     month + DISTID | 0 | DISTID, data = panel_af)
IQ.m.g.t.q <- felm(unconstrained1 ~ temperature + temperature_squared + precipitation + constrained1 + ied +
                     temperature_1 + temperature_2 + temperature_3 + temperature_4 + temperature_5 + temperature_6 + temperature_7 +
                     precipitation_1 + precipitation_2 + precipitation_3 + precipitation_4 + precipitation_5 + precipitation_6 + precipitation_7 +
                     unconstrained1_1 + unconstrained1_2 + unconstrained1_3 + unconstrained1_4 + unconstrained1_5 + unconstrained1_6 + unconstrained1_7 +
                     constrained1_1 + constrained1_2 + constrained1_3 + constrained1_4 + constrained1_5 + constrained1_6 + constrained1_7 +
                     ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                     month + ADM3NAME | 0 | ADM3NAME, data = panel_iq)

AF.m.g.t.c <- felm(df ~ temperature + temperature_squared + temperature_cubed + precipitation + idf + ied +
                     temperature_1 + temperature_2 + temperature_3 + temperature_4 + temperature_5 + temperature_6 + temperature_7 +
                     precipitation_1 + precipitation_2 + precipitation_3 + precipitation_4 + precipitation_5 + precipitation_6 + precipitation_7 +
                     df_1 + df_2 + df_3 + df_4 + df_5 + df_6 + df_7 +
                     idf_1 + idf_2 + idf_3 + idf_4 + idf_5 + idf_6 + idf_7 +
                     ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                     month + DISTID | 0 | DISTID, data = panel_af)
IQ.m.g.t.c <- felm(unconstrained1 ~ temperature + temperature_squared + temperature_cubed + precipitation + constrained1 + ied +
                     temperature_1 + temperature_2 + temperature_3 + temperature_4 + temperature_5 + temperature_6 + temperature_7 +
                     precipitation_1 + precipitation_2 + precipitation_3 + precipitation_4 + precipitation_5 + precipitation_6 + precipitation_7 +
                     unconstrained1_1 + unconstrained1_2 + unconstrained1_3 + unconstrained1_4 + unconstrained1_5 + unconstrained1_6 + unconstrained1_7 +
                     constrained1_1 + constrained1_2 + constrained1_3 + constrained1_4 + constrained1_5 + constrained1_6 + constrained1_7 +
                     ied_1 + ied_2 + ied_3 + ied_4 + ied_5 + ied_6 + ied_7 |
                     month + ADM3NAME | 0 | ADM3NAME, data = panel_iq)

###
### Create dataset for prediction:
###

new_data_af <- data.frame(temperature = temperature_range, 
                          temperature_squared = (temperature_range)^2,
                          temperature_cubed = (temperature_range)^3,
                          
                          precipitation = mean(panel_af$precipitation, na.rm = T), 
                          df = mean(panel_af$df), 
                          idf = mean(panel_af$idf),  
                          ied = mean(panel_af$ied), 
                          
                          temperature_1 = mean(panel_af$temperature_1, na.rm = T),
                          temperature_2 = mean(panel_af$temperature_2, na.rm = T),
                          temperature_3 = mean(panel_af$temperature_3, na.rm = T),
                          temperature_4 = mean(panel_af$temperature_4, na.rm = T),
                          temperature_5 = mean(panel_af$temperature_5, na.rm = T),
                          temperature_6 = mean(panel_af$temperature_6, na.rm = T),
                          temperature_7 = mean(panel_af$temperature_7, na.rm = T),
                          
                          precipitation_1 = mean(panel_af$precipitation_1, na.rm = T),
                          precipitation_2 = mean(panel_af$precipitation_2, na.rm = T),
                          precipitation_3 = mean(panel_af$precipitation_3, na.rm = T),
                          precipitation_4 = mean(panel_af$precipitation_4, na.rm = T),
                          precipitation_5 = mean(panel_af$precipitation_5, na.rm = T),
                          precipitation_6 = mean(panel_af$precipitation_6, na.rm = T),
                          precipitation_7 = mean(panel_af$precipitation_7, na.rm = T),
                          
                          df_1 = mean(panel_af$df_1, na.rm = T),
                          df_2 = mean(panel_af$df_2, na.rm = T),
                          df_3 = mean(panel_af$df_3, na.rm = T),
                          df_4 = mean(panel_af$df_4, na.rm = T),
                          df_5 = mean(panel_af$df_5, na.rm = T),
                          df_6 = mean(panel_af$df_6, na.rm = T),
                          df_7 = mean(panel_af$df_7, na.rm = T),
                          
                          idf_1 = mean(panel_af$idf_1, na.rm = T),
                          idf_2 = mean(panel_af$idf_2, na.rm = T),
                          idf_3 = mean(panel_af$idf_3, na.rm = T),
                          idf_4 = mean(panel_af$idf_4, na.rm = T),
                          idf_5 = mean(panel_af$idf_5, na.rm = T),
                          idf_6 = mean(panel_af$idf_6, na.rm = T),
                          idf_7 = mean(panel_af$idf_7, na.rm = T),
                          
                          ied_1 = mean(panel_af$ied_1, na.rm = T),
                          ied_2 = mean(panel_af$ied_2, na.rm = T),
                          ied_3 = mean(panel_af$ied_3, na.rm = T),
                          ied_4 = mean(panel_af$ied_4, na.rm = T),
                          ied_5 = mean(panel_af$ied_5, na.rm = T),
                          ied_6 = mean(panel_af$ied_6, na.rm = T),
                          ied_7 = mean(panel_af$ied_7, na.rm = T),
                          
                          week = median(panel_af$week),
                          month = median(panel_af$month),
                          DISTID = 1001)

##

new_data_iq <- data.frame(temperature = temperature_range, 
                          temperature_squared = (temperature_range)^2,
                          temperature_cubed = (temperature_range)^3,
                          
                          precipitation = mean(panel_iq$precipitation, na.rm = T), 
                          unconstrained1 = mean(panel_iq$unconstrained1), 
                          constrained1 = mean(panel_iq$constrained1),  
                          ied = mean(panel_iq$ied),  
                          
                          temperature_1 = mean(panel_iq$temperature_1, na.rm = T),
                          temperature_2 = mean(panel_iq$temperature_2, na.rm = T),
                          temperature_3 = mean(panel_iq$temperature_3, na.rm = T),
                          temperature_4 = mean(panel_iq$temperature_4, na.rm = T),
                          temperature_5 = mean(panel_iq$temperature_5, na.rm = T),
                          temperature_6 = mean(panel_iq$temperature_6, na.rm = T),
                          temperature_7 = mean(panel_iq$temperature_7, na.rm = T),
                          
                          precipitation_1 = mean(panel_iq$precipitation_1, na.rm = T),
                          precipitation_2 = mean(panel_iq$precipitation_2, na.rm = T),
                          precipitation_3 = mean(panel_iq$precipitation_3, na.rm = T),
                          precipitation_4 = mean(panel_iq$precipitation_4, na.rm = T),
                          precipitation_5 = mean(panel_iq$precipitation_5, na.rm = T),
                          precipitation_6 = mean(panel_iq$precipitation_6, na.rm = T),
                          precipitation_7 = mean(panel_iq$precipitation_7, na.rm = T),
                          
                          unconstrained1_1 = mean(panel_iq$unconstrained1_1, na.rm = T),
                          unconstrained1_2 = mean(panel_iq$unconstrained1_2, na.rm = T),
                          unconstrained1_3 = mean(panel_iq$unconstrained1_3, na.rm = T),
                          unconstrained1_4 = mean(panel_iq$unconstrained1_4, na.rm = T),
                          unconstrained1_5 = mean(panel_iq$unconstrained1_5, na.rm = T),
                          unconstrained1_6 = mean(panel_iq$unconstrained1_6, na.rm = T),
                          unconstrained1_7 = mean(panel_iq$unconstrained1_7, na.rm = T),
                          
                          constrained1_1 = mean(panel_iq$constrained1_1, na.rm = T),
                          constrained1_2 = mean(panel_iq$constrained1_2, na.rm = T),
                          constrained1_3 = mean(panel_iq$constrained1_3, na.rm = T),
                          constrained1_4 = mean(panel_iq$constrained1_4, na.rm = T),
                          constrained1_5 = mean(panel_iq$constrained1_5, na.rm = T),
                          constrained1_6 = mean(panel_iq$constrained1_6, na.rm = T),
                          constrained1_7 = mean(panel_iq$constrained1_7, na.rm = T),
                          
                          ied_1 = mean(panel_iq$ied_1, na.rm = T),
                          ied_2 = mean(panel_iq$ied_2, na.rm = T),
                          ied_3 = mean(panel_iq$ied_3, na.rm = T),
                          ied_4 = mean(panel_iq$ied_4, na.rm = T),
                          ied_5 = mean(panel_iq$ied_5, na.rm = T),
                          ied_6 = mean(panel_iq$ied_6, na.rm = T),
                          ied_7 = mean(panel_iq$ied_7, na.rm = T),
                          
                          week = median(panel_iq$week),
                          month = median(panel_iq$month),
                          ADM3NAME = "Zakho")

###
### Then, get predictions:
###

# The following code is sourced from Stack Overflow:
# https://stackoverflow.com/questions/30491545/predict-method-for-felm-from-lfe-package

predict.felm <- function(object, newdata, se.fit = FALSE,
                         interval = "none",
                         level = 0.95){
  if(missing(newdata)){
    stop("predict.felm requires newdata and predicts for all group effects = 0.")
  }
  
  tt <- terms(object)
  Terms <- delete.response(tt)
  attr(Terms, "intercept") <- 0
  
  m.mat <- model.matrix(Terms, data = newdata)
  m.coef <- as.numeric(object$coef)
  fit <- as.vector(m.mat %*% object$coef)
  fit <- data.frame(fit = fit)
  
  if(se.fit | interval != "none"){
    if(!is.null(object$clustervcv)){
      vcov_mat <- object$clustervcv
    } else if (!is.null(object$robustvcv)) {
      vcov_mat <- object$robustvcv
    } else if (!is.null(object$vcv)){
      vcov_mat <- object$vcv
    } else {
      stop("No vcv attached to felm object.")
    }
    se.fit_mat <- sqrt(diag(m.mat %*% vcov_mat %*% t(m.mat)))
  }
  if(interval == "confidence"){
    t_val <- qt((1 - level) / 2 + level, df = object$df.residual)
    fit$lwr <- fit$fit - t_val * se.fit_mat
    fit$upr <- fit$fit + t_val * se.fit_mat
  } else if (interval == "prediction"){
    stop("interval = \"prediction\" not yet implemented")
  }
  if(se.fit){
    return(list(fit=fit, se.fit=se.fit_mat))
  } else {
    return(fit)
  }
}

# Afghanistan:
pred_AF.w.g.t.l <- predict.felm(AF.w.g.t.l, new_data_af)
pred_AF.w.g.t.q <- predict.felm(AF.w.g.t.q, new_data_af)
pred_AF.w.g.t.c <- predict.felm(AF.w.g.t.c, new_data_af)

pred_AF.m.g.t.l <- predict.felm(AF.m.g.t.l, new_data_af)
pred_AF.m.g.t.q <- predict.felm(AF.m.g.t.q, new_data_af)
pred_AF.m.g.t.c <- predict.felm(AF.m.g.t.c, new_data_af)

# Iraq:
pred_IQ.w.g.t.l <- predict.felm(IQ.w.g.t.l, new_data_iq)
pred_IQ.w.g.t.q <- predict.felm(IQ.w.g.t.q, new_data_iq)
pred_IQ.w.g.t.c <- predict.felm(IQ.w.g.t.c, new_data_iq)

pred_IQ.m.g.t.l <- predict.felm(IQ.m.g.t.l, new_data_iq)
pred_IQ.m.g.t.q <- predict.felm(IQ.m.g.t.q, new_data_iq)
pred_IQ.m.g.t.c <- predict.felm(IQ.m.g.t.c, new_data_iq)

###
### Then, build dataset from the predictions (for statistically significant models only:
###

preds_ols <- cbind(pred_AF.w.g.t.l, pred_AF.w.g.t.c, pred_AF.m.g.t.l, pred_AF.m.g.t.c,
                   pred_IQ.w.g.t.q, pred_IQ.w.g.t.c, pred_IQ.m.g.t.q)
colnames(preds_ols) <- c("model1_af1", "model2_af2", "model3_af3", "model4_af4",
                         "model5_iq1", "model6_iq2", "model7_iq3")
# Create temperature range:
preds_ols$temperature <- temperature_range

# Rescale for purposes of comparability:

preds_ols[, 1:4] <- preds_ols[, 1:4] * 20
preds_ols[, 5:7] <- preds_ols[, 5:7] * 6

preds_ols$mean_af <- rowMeans(preds_ols[, 1:4])
preds_ols$mean_iq <- rowMeans(preds_ols[, 5:7])

# Then, distinguish predictions for temperatures outside of each country's observed values:
temperature_range_af <- round(min(panel_af$temperature)):round(max(panel_af$temperature))
temperature_range_iq <- round(min(panel_iq$temperature)):round(max(panel_iq$temperature))

preds_ols$model1_af1_c <- preds_ols$model2_af2_c <- preds_ols$model3_af3_c <- preds_ols$model4_af4_c <-
  preds_ols$model5_iq1_c <- preds_ols$model6_iq2_c <- preds_ols$model7_iq3_c <- 
  preds_ols$mean_af_c <- preds_ols$mean_iq_c <- NA

preds_ols[which(preds_ols$temperature %in% temperature_range_af), ]$model1_af1_c <- preds_ols[which(preds_ols$temperature %in% temperature_range_af), ]$model1_af1
preds_ols[which(preds_ols$temperature %in% temperature_range_af), ]$model2_af2_c <- preds_ols[which(preds_ols$temperature %in% temperature_range_af), ]$model2_af2
preds_ols[which(preds_ols$temperature %in% temperature_range_af), ]$model3_af3_c <- preds_ols[which(preds_ols$temperature %in% temperature_range_af), ]$model3_af3
preds_ols[which(preds_ols$temperature %in% temperature_range_af), ]$model4_af4_c <- preds_ols[which(preds_ols$temperature %in% temperature_range_af), ]$model4_af4
preds_ols[which(preds_ols$temperature %in% temperature_range_iq), ]$model5_iq1_c <- preds_ols[which(preds_ols$temperature %in% temperature_range_iq), ]$model5_iq1
preds_ols[which(preds_ols$temperature %in% temperature_range_iq), ]$model6_iq2_c <- preds_ols[which(preds_ols$temperature %in% temperature_range_iq), ]$model6_iq2
preds_ols[which(preds_ols$temperature %in% temperature_range_iq), ]$model7_iq3_c <- preds_ols[which(preds_ols$temperature %in% temperature_range_iq), ]$model7_iq3

preds_ols[which(preds_ols$temperature %in% temperature_range_af), ]$mean_af_c <- preds_ols[which(preds_ols$temperature %in% temperature_range_af), ]$mean_af
preds_ols[which(preds_ols$temperature %in% temperature_range_iq), ]$mean_iq_c <- preds_ols[which(preds_ols$temperature %in% temperature_range_iq), ]$mean_iq

###
### Then, plot:
###

p1 <- ggplot() + 
  geom_histogram(data = all_observed_temps[all_observed_temps$temperature %in% temperature_range, , drop=F], aes(x=temperature, y=(..ncount..)), bins = 100, size=.1, color="black", alpha=.4) +
  
  geom_line(data = preds_ols, aes(x = temperature, y = model1_af1), color = "blue4", size = 0.15, linetype = "dashed") +
  geom_line(data = preds_ols, aes(x = temperature, y = model2_af2), color = "blue4", size = 0.15, linetype = "dotted") +
  geom_line(data = preds_ols, aes(x = temperature, y = model3_af3), color = "blue4", size = 0.15, linetype = "dashed") +
  geom_line(data = preds_ols, aes(x = temperature, y = model4_af4), color = "blue4", size = 0.15, linetype = "dotted") +
  geom_line(data = preds_ols, aes(x = temperature, y = model5_iq1), color = "darkred", size = 0.15, linetype = "dashed") +
  geom_line(data = preds_ols, aes(x = temperature, y = model6_iq2), color = "darkred", size = 0.15, linetype = "dotted") +
  geom_line(data = preds_ols, aes(x = temperature, y = model7_iq3), color = "darkred", size = 0.15, linetype = "dashed") +

  geom_line(data = preds_ols, aes(x = temperature, y = model1_af1_c), color = "blue4", size = 0.5, linetype = "dashed") +
  geom_line(data = preds_ols, aes(x = temperature, y = model2_af2_c), color = "blue4", size = 0.5, linetype = "dotted") +
  geom_line(data = preds_ols, aes(x = temperature, y = model3_af3_c), color = "blue4", size = 0.5, linetype = "dashed") +
  geom_line(data = preds_ols, aes(x = temperature, y = model4_af4_c), color = "blue4", size = 0.5, linetype = "dotted") +
  geom_line(data = preds_ols, aes(x = temperature, y = model5_iq1_c), color = "darkred", size = 0.5, linetype = "dashed") +
  geom_line(data = preds_ols, aes(x = temperature, y = model6_iq2_c), color = "darkred", size = 0.5, linetype = "dotted") +
  geom_line(data = preds_ols, aes(x = temperature, y = model7_iq3_c), color = "darkred", size = 0.5, linetype = "dashed") +

  # geom_line(data = preds_ols, aes(x = temperature, y = mean), color = "black", size = 2) +
  geom_line(data = preds_ols, aes(x = temperature, y = mean_af_c), color = "blue4", size = 1) +
  geom_line(data = preds_ols, aes(x = temperature, y = mean_iq_c), color = "darkred", size = 1) +

  geom_line(data = preds_ols, aes(x = temperature, y = mean_af), color = "blue4", size = 0.5) +
  geom_line(data = preds_ols, aes(x = temperature, y = mean_iq), color = "darkred", size = 0.5) +

  theme_bw() +
  geom_vline(xintercept = preds_ols[which(preds_ols$mean_af == max(preds_ols$mean_af)), ]$temperature-0.25, linetype=4, color = "blue4") +
  geom_vline(xintercept = preds_ols[which(preds_ols$mean_iq == max(preds_ols$mean_iq)), ]$temperature, linetype=4, color = "darkred") +

  ggtitle("") +
  ylab("Predicted Instances of Violent Attacks") + 
  xlab("Temperature") +
  xlim(-20,120) + 
  ylim(0,2.5) + 
  
  scale_x_continuous(breaks = seq(-20, 120, by = 5)) + 
  theme(axis.text.y = element_blank(), 
        axis.ticks = element_blank(), 
        text = element_text(size=10),
        axis.text.x = element_text(angle = 45)) +
  theme(axis.text.x = element_text(face="bold", size=12),
        axis.text.y = element_text(face="bold",
                                   size=14),
        axis.title=element_text(size=14,face="bold"))

# PRODUCES FIGURE 10 IN PAPER:
p1

###
### Then, report the results in a table:
###

stargazer(AF.w.g.t.l, AF.w.g.t.c, AF.m.g.t.l, AF.m.g.t.c,
          IQ.w.g.t.q, IQ.w.g.t.c, IQ.m.g.t.q, 
          model.names = TRUE, 
          omit.yes.no = c("Yes", "No"), 
          star.cutoffs = c(.1, .05, .01),
          omit = "factor", single.row = F)

###
### END.
###
